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Abstract. Deterministic chaotic dynamics presumes that the state space can be partitioned arbi- 
trarily finely. In a physical system, the inevitable presence of some noise sets a finite limit to the 
finest possible resolution that can be attained. Much previous research deals with what this attain- 
able resolution might be, all of it based on a global averages over stochastic flow. We show how 
to compute the locally optimal partition, for a given dynamical system and given noise, in terms of 
local eigenfunctions of the Fokker-Planck operator and its adjoint. We first analyze the interplay of 
the deterministic dynamics with the noise in the neighborhood of a periodic orbit of a map, by us- 
ing a discretized version of Fokker-Planck formalism. Then we propose a method to determine the 
'optimal resolution' of the state space, based on solving Fokker-Planck's equation locally, on sets 
of unstable periodic orbits of the deterministic system. We test our hypothesis on unimodal maps. 
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1. INTRODUCTION 

The effect of noise on the behavior of a nonlinear dynamical system is a fundamental 
problem in many areas of science [1, 2, 3], and the interplay of noise and chaotic 
dynamics is of particular interest [4, 5, 6]. Our purpose here is two-fold. First, we address 
operationally the fact that weak noise limits the attainable resolution of the state space 
of a chaotic system by formulating the optimal partition hypothesis. In ref. [7] we have 
shown that the hypothesis enables us to define the optimal partition for a 1 -dimensional 
map; here we explain how it is implemented for a high-dimensional state space flows 
with a few expanding directions, such as the transitional Re number Navier-Stokes flows. 
Second, we show that the optimal partition hypothesis replaces the Fokker-Planck PDEs 
by finite, low-dimensional matrix Fokker-Planck operators, with finite cycle expansions, 
optimal for a given level of precision, and whose eigenvalues give good estimates of 
long-time observables (escape rates, Lyapunov exponents, etc.). 

A chaotic trajectory explores a strange attractor, and for chaotic flows evaluation 
of long-time averages requires effective partitioning of the state space into smaller re- 
gions. In a hyperbolic, everywhere unstable deterministic dynamical system, consecutive 
Poincare section returns subdivide the state space into exponentially growing number of 
regions, each region labeled by a distinct finite symbol sequence, as in figure 1. In the 
unstable directions these regions stretch, while in the stable directions they shrink ex- 
ponentially. The set of unstable periodic orbits forms a 'skeleton' that can be used to 
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FIGURE 1. (a) A coarse partition of state space ^ into regions «^o> M\, and ^#2, labeled by ternary 
alphabet srf = {1,2,3}. (b) A 1-step memory refinement of the partition of figure 1, with each region ^ 
subdivided into > an d *^2> labeled by nine 'words' {00, 01 , 02, • • • , 21 , 22}. 

implement such partition of the state space, each region a neighborhood of a periodic 
point [8, 9]. Longer and longer cycles yield finer and finer partitions as the neighbor- 
hood of each unstable cycle p shrinks exponentially with cycle period as l/|A p |, where 
A p is the product of cycle's expanding Floquet multipliers. As there is an exponen- 
tially growing infinity of longer and longer cycles, with each neighborhood shrinking 
asymptotically to a point, a deterministic chaotic system can - in principle - be resolved 
arbitrarily finely. But that is a fiction for any of the following reasons: 

• any physical system experiences (background, observational, intrinsic, measure- 
ment, ■ ■ ■ ) noise 

• any numerical computation is a noisy process due to the finite precision of each 
step of computation 

• any set of dynamical equations models nature up to a given finite accuracy, since 
degrees of freedom are always neglected 

• any prediction only needs to be computed to a desired finite accuracy 

The problem we address here is sketched in figure 2; while a deterministic partition 
can, in principle, be made arbitrarily fine, in practice any noise will blur the boundaries 
and render the best possible partition finite. Thus our task is to determine the optimal 
attainable resolution of the state space of a given hyperbolic dynamical system, affected 
by a given weak noise. This we do by formulating the optimal partition hypothesis 
which we believe determines the best possible state space partition for a desired level of 
predictive precision. We know of no practical way of computing the 'blurred' partition 
boundaries of figure 2(b). Instead, we propose to determine the optimal partition in 
terms of blurring of periodic point neighborhoods, as in figure 2 (d). As we demonstrate 
in sect. 5, our implementation requires determination of only a small set of solutions of 
the deterministic equations of motion. 

Intuitively, the noise smears out the neighborhood of a periodic point, whose size is 
now determined by the interplay between the diffusive spreading parameterized [10, 11, 
12] by a diffusion constant, and its exponentially shrinking deterministic neighborhood. 
If the noise is weak, the short-time dynamics is not altered significantly: short periodic 
orbits of the deterministic flow still coarsely partition the state space. As the periods 
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FIGURE 2. (a) A deterministic partition of state space jtft '. (b) Noise blurs the partition boundaries. At 
some level of deterministic partitioning some of the boundaries start to overlap, preventing further local 
refinement of adjecent neighborhoods, (c) The fixed point 1 = {x{\ and the two-cycle 01 = {jtoii-^io} are 
examples of the shortest period periodic points within the partition regions of (a). The optimal partition 
hypothesis: (d) Noise blurs periodic points into cigar-shaped trajectory -centered densities explored by the 
Langevin noise. The optimal partition hypothesized in this paper consists of the maximal set of resolvable 
periodic point neighborhoods. 



of periodic orbits increase, the diffusion always wins, and successive refinements of a 
deterministic partition of the state space stop at the finest attainable partition, beyond 
which the diffusive smearing exceeds the size of any deterministic subpartition. 

There is a considerable literature (reviewed here in remark 5.1) on interplay of noise 
and chaotic deterministic dynamics, and the closely related problem of limits on the va- 
lidity of the semi-classical periodic orbit quantization. All of this literature implicitly 
assumes uniform hyperbolicity and seeks to define a single, globally averaged, diffusion 
induced average resolution (Heisenberg time, in the context of semi-classical quanti- 
zation). However, the local diffusion rate differs from a trajectory to a trajectory, as 
different neighborhoods merge at different times, so there is no one single time beyond 
which noise takes over. Nonlinear dynamics interacts with noise in a nonlinear way, and 
methods for implementing the optimal partition for a given noise still need to be devel- 
oped. This paper is an attempt in this direction. Here we follow and expand upon the 
Fokker-Planck approach to the 'optimal partition hypothesis' introduced in ref. [7]. 

What is novel here is that we show how to compute the locally optimal partition, 
for a given dynamical system and given noise, in terms of local eigenfunctions of the 
forward-backward actions of the Fokker-Planck operator and its adjoint. This is much 
simpler than it sounds: the Lyapunov equation 



Q = MQM T + A 



(and its generalizations to periodic points of hyperbolic flows), determines Q, the size 
of local neighborhood, as balance of the noise variance A and the linearized dynamics 
M. The effort of going local brings a handsome reward: as the optimal partition is 
always finite, the dynamics on this 'best possible of all partitions' is encoded by a finite 
transition graph of finite memory, and the Fokker- Planck operator can be represented by 
a finite matrix. In addition, while the state space of a generic deterministic flow is an 
infinitely interwoven hierarchy of attracting, hyperbolic, elliptic and parabolic regions, 
the noisy dynamics erases any structures finer than the optimal partition, thus curing both 
the affliction of long-period attractors/elliptic islands with very small immediate basins 
of attraction/ellipticity, and the power-law correlation decays caused by marginally 
stable regions of state space. 

The dynamical properties of high-dimensional flows are not just simple extensions 
of lower-dimensional dynamics, and a persuasive application of the Ruelle / Gutzwiller 
periodic orbit theory to high-dimensional dynamics would be an important advance. If 
such flow has only a few expanding directions, the above set of overlapping stochastic 
'cigars' should provide an optimal, computable cover of the long-time chaotic attractor 
embedded in a state space of arbitrarily high dimension. 

The requisite Langevin / Fokker-Planck description of noisy flows is reviewed in 
sect. 2. This discussion leans heavily on the deterministic dynamics and periodic orbit 
theory notation, summarized in appendix A. In sect. 3 we derive the formulas for the 
size of noise-induced neighborhoods of attractive fixed and periodic points, for maps 
and flows in arbitrary dimension. These formulas are known as Lyapunov equations, 
reviewed in appendix C. In order to understand the effect on noise on the hyperbolic, 
mixed expanding / contracting dynamics, we study the eigenfunctions of the discrete 
time Fokker-Planck operator in linear neighborhood of a fixed point of a noisy one- 
dimensional map in sect. 4, and show that the neighborhood along unstable directions is 
fixed by the evolution of a Gaussian density of trajectories under the action of the adjoint 
Fokker-Planck operator. The continuous time formulation of the same problem, known 
as the Ornstein-Uhlenbeck process, is reviewed in appendix B. Having defined the local 
neighborhood of every periodic point, we turn to the global partition problem. Previous 
attempts at state space partitioning are reviewed in remark 5.1. We formulate our optimal 
partition hypothesis in sect. 5: track the diffusive widths of unstable periodic orbits until 
they start to overlap. We test the approach by applying it to a 1 -dimensional repeller, 
and in sect. 6, we assess the accuracy of our method by computing the escape rate and 
the Lyapunov exponent, discuss weak noise corrections, and compare the results with a 
discretization of the Fokker-Planck operator on a uniform mesh. In sect. 7 we address the 
problem of estimating the optimal partition of a non-hyperbolic map, where the linear 
approximation to the Fokker-Planck operator fails. The results are summarized and the 
open problems discussed in sect. 8. 

2. NOISY TRAJECTORIES AND THEIR DENSITIES 

The literature on stochastic dynamical systems is vast, starting with the Laplace 1810 
memoir [13]. The material reviewed in this section, sect. 4 and appendix B is standard [1, 
3, 14], but needed in order to set the notation for what is new here, the role that Fokker- 



Planck operators play in defining stochastic neighborhoods of periodic orbits. The key 
result derived here is the well known evolution law (14) for the covariance matrix Q a of 
a linearly evolved Gaussian density, 

Qa+i =M a Q a M T a +K. 

To keep things simple we shall use only the discrete time dynamics in what follows, but 
we do discuss the continuous time formulation in appendix B, as our results apply both 
to the continuous and discrete time flows. 

Consider a noisy discrete time dynamical system [15, 16, 17] 

X n +1 = f{x n ) + , (1) 

where x is a J-dimensional state vector, and x n j is its jth component at time n. In the 
Fokker-Planck description individual noisy trajectories are replaced by the evolution of 
the density of noisy trajectories, with the x n+ i — f(x n ) probability distribution of zero 
mean and covariance matrix (diffusion tensor) A, 

(Snj)=0, =A, 7 S nm , (2) 

where (■ ■ ■) stands for ensemble average over many realizations of the noise. 

The general case of a diffusion tensor A(x) which is a state space position dependent 
but time independent can be treated along the same lines. In this case the stochastic flow 
(1) is written as x n+ \ = x n + o{x) t, n , where (£ n ^} = 1 8 nm is white noise, A = o o T , 
o{x) is called the 'diffusion matrix', and the noise is referred to as 'multiplicative' (see 
Kuehn[18]). 

The action of discrete one-time step Fokker-Planck operator on the density distribu- 
tion p at time k, 

Pk+l(y) = l^Pk](y) = J dx^(y,x)p k (x) 

J?(y,x) = Ie4(>'-/W) r i(y-/W), (3 ) 

is centered on the deterministic step f(x) and smeared out diffusively by noise. Were 
diffusion uniform and isotropic, A(^) = 2D1, the Fokker-Planck operator would be pro- 
portional to exp (— {y — f(x)} 2 /2A) , i.e., the penalty for straying from the deterministic 
path is just a quadratic error function. The kth iterate of is a J-dimensional path 
integral over the k — 1 intermediate noisy trajectory points, 

Sf k (x k ,xo) = J [dx] e -^ n+ i-fM) T ^n +l -fM) ^ (4) 

where the Gaussian normalization factor in (3) is absorbed into intermediate integrations 
by defining 

k-l j d 

[ J *] = n^> N= v2^detA. (5) 



We shall also need to determine the effect of noise accumulated along the trajectory 
points preceding x. As the noise is additive forward in time, one cannot simply invert 
the Fokker-Planck operator; instead, the past is described by the adjoint Fokker-Planck 
operator, 



which transports a density concentrated around the point f(x) to a density concentrated 
around the previous point x and adds noise to it. In the deterministic, vanishing noise 
limit this is the Koopman operator (76). 

The Fokker-Planck operator (3) is non-hermitian and non-unitary. For example, if the 
deterministic flow is contracting, the natural measure (the leading right eigenvector of 
the Fokker-Planck operator) will be concentrated and peaked, but then the corresponding 
left eigenvector has to be broad and flat, as backward in time the deterministic flow 
is expanding. We shall denote by p a the right eigenvectors of and by p a its left 
eigenvectors, i.e., the right eigenvectors of the adjoint operator ' . 



Our first goal is to convince the reader that the diffusive dynamics of nonlinear flows is 
fundamentally different from. Brownian motion, with the flow inducing a local, history 
dependent noise. In order to accomplish this, we go beyond the standard stochastic 
literature and generalize the notion of invariant deterministic recurrent solutions, such as 
fixed points and periodic orbits, to noisy flows. While a Langevin trajectory (1) cannot 
be periodic, in the Fokker-Planck formulation (4) a recurrent motion can be defined 
as one where a peaked distribution returns to the initial neighborhood after time n. 
Recurrence so defined not only coincides with the classical notion of a recurrent orbit 
in the vanishing noise limit, but it also enables us to derive exact formulas for how this 
local, history dependent noise is to be computed. 

As the function x n+ \ — f(x n ) is a nonlinear function, in general the path integral (4) 
can only be evaluated numerically. In the vanishing noise limit the Gaussian kernel 
sharpens into the Dirac 8 -function, and the Fokker-Planck operator reduces to the 
deterministic Perron-Frobenius operator (75). For weak noise the Fokker-Planck oper- 
ator can be evaluated perturbatively [19, 20, 21, 22] as an asymptotic series in powers of 
the diffusion constant, centered on the deterministic trajectory. Here we retain only the 
linear term in this series, which has a particulary simple dynamics given by a covariance 
matrix evolution formula (see (14) and (22) below) that we now derive. 

We shift local coordinates to the deterministic trajectory xq, x\, X2,---,} 

centered coordinate frames x = x a +z a , Taylor expand f(x) = f a (z a ) = x a +i +M a z a H , 

and approximate the noisy map (1) by its linearization, 



with the deterministic trajectory points at z a = z a +i = 0, and M a the one step Jacobian 
matrix (see appendix A). The corresponding linearized Fokker-Planck operator (3) is 




(6) 



3. ALL NONLINEAR NOISE IS LOCAL 



Za+l = M aZa + %a , M i} = df/dxj 



(7) 



given in the local coordinates p a (z a ) = p(x a + z a ,a) by 

Pa+l(z a +l) = J dz a ^ a{Za+\,Z a ) Pa{Za) (8) 

by the linearization (7) centered on the deterministic trajectory 

&a{Za+\,Za) = ig-^-Wl (^i-M*,) . (9) 

The subscript 'a' in a distinguishes the local, linearized Fokker-Planck operator from 
the full operator (4). 

The kernel of the linearized Fokker-Planck operator (9) is a Gaussian. As a convo- 
lution of a Gaussian with a Gaussian is again a Gaussian, we investigate the action of 
the linearized Fokker-Planck operator on a normalized, cigar-shaped Gaussian density 
distribution 

Pa{z) = ^e-^Wa\ Ca = (2K) d l 2 (tetQ a ) l l\ (10) 
and the action of the linearized adjoint Fokker-Planck operator on density 

Pa{z) = ^-«"^± Z , C a = (2K) d / 2 (daQ a ) l l\ (11) 

also centered on the deterministic trajectory, but with its own strictly positive [dxd] 
covariance matrices Q, Q. Label 'a' plays a double role, and {a + I, a} stands both for 
the {next, initial} space partition and for the times the trajectory lands in these partitions 
(see appendix A. 1). The linearized Fokker-Planck operator (9) maps the Gaussian p a (z a ) 
into the Gaussian 

Pa+l(z a+ l) = 1 (12) 

one time step later. Likewise, linearizing the adjoint Fokker-Planck operator (6) around 
the x a trajectory point yields: 

. , , 1 fr , , -\[(Za+\-M a Za) T ^{z a+ l-M a Za)+Z T a+l -J— Z a +l] 

Completing the squares, integrating and substituting (10), respectively (11) we obtain 
the formula for Q covariance matrix evolution forward in time, 

Q a+ i=M a Q a M T a +K, (14) 

and in the adjoint case, the evolution of the Q is given by 

M a Q a M T a =Q a+l +K. (15) 

The two covariance matrices differ, as the adjoint evolution Q a is computed by going 
backwards along the trajectory. These covariance evolution rules are the basis of all that 
follows, except for the 'flat-top' of sect. 7. 



Think of the initial co variance matrix (10) as an error matrix describing the precision 
of the initial state, a cigar-shaped probability distribution p a (z a )- In one time step this 
density is deterministically advected and deformed into density with covariance MQM T , 
and then the noise A is added: the two kinds of independent uncertainties add up as 
sums of squares, hence the covariance evolution law (14), resulting in the Gaussian 
ellipsoid whose widths and orientation are given by the singular values and singular 
vectors (68) of the covariance matrix. After n time steps, the variance Q a is built up 
from the deterministically propagated M"Q a - n M" T initial distribution, and the sum of 
noise kicks at intervening times, M^A^^M^ 7 , also propagated deterministically. 

The pleasant surprise is that the evaluation of this noise requires no Fokker-Planck 
PDE formalism. The width of a Gaussian packet centered on a trajectory is fully spec- 
ified by a deterministic computation that is already a pre-computed byproduct of the 
periodic orbit computations; the deterministic orbit and its linear stability. We have at- 
tached label 'a' to A a = A(x a ) in (14) to account for the noise distributions that are 
inhomogeneous, state space dependent, but time independent multiplicative noise. As 
we shall show, in nonlinear dynamics the noise is never isotropic and/or homogeneous. 
For example, if the iterative system we are studying is obtained by Poincare sections of a 
continuous time flow, the accumulated noise integrated over one Poincare section return 
depends on the return trajectory segment, even when the infinitesimal time step noise 
(87) is homogenous. 

Remark 3.1 Covariance evolution. In quantum mechanics the linearized evolution operator 
corresponding to the linearized Fokker-Planck operator (9) is known as the Van Vleck propa- 
gator, the basic block in the semi-classical periodic orbit quantization [23, 24]. Q covariance 
matrix composition rule (14) or its continuous time version (106) is called 'covariance evolu- 
tion' in ref. [25], for example, but it goes all the way back to Lyapunov's 1892 thesis [26], see 
appendix C. In the Kalman filter literature [27, 28] it is called 'prediction'. 



3.1. The attractive, the repulsive and the noisy 

For Browniam dynamics x n+ i = x n + £ n , with M = 1, we obtain Q n = Qq + nA, i.e., 
the variance of a Gaussian packet of p n (z) noisy trajectories grows linearly in time, as 
expected for the Brownian diffusion. What happens for nontrivial, M/l, dynamics? 
The formulas (14) and (15) are exact for finite numbers of time steps, but whether they 
have a long time limit depends on the stability of the deterministic trajectory. 

Here we shall derive the n — > oo limit for deterministic flows that are either contracting 
or expanding in all eigen-directions, with asymptotic stationary distributions concen- 
trated either on fixed points or periodic points. We shall consider the general hyperbolic 
flows, where some of the eigen-directions are unstable, and other stable in another pub- 
lication [29]. In this context the description in terms of periodic orbits is very useful; the 
neighborhood of a periodic point will be defined as the noise contracting neighborhood 
forward in time along contracting eigen-directions, backward in time along the unstable, 
expanding eigen-directions. The short cycles will be the most important ones, and only 
finite time, single cycle period calculations will be required. 



If M is contracting, with the multipliers 1 > |Ai | > | A.2I > . . . > |A^| , in n time steps 
the memory of the covariance Q a - n of the starting density is forgotten at exponential 
rate ~ |Aj |~ 2 ", with iteration of (14) leading to a limit distribution: 

Q a = A a +M a _ 1 A a _ 1 Mj_ 1 + M 2 _ 2 A fl _ 2 (M 2 _ 2 ) r + . . . . (16) 

For fixed and periodic points we can give an explicit formula for the n — > 00 covariance. 

Consider a noisy map (1) with a deterministic fixed point at x q . In a neighborhood 
x = x q + z we approximate the map / by its linearization (7) with the fixed point at z = 0, 
acting on a Gaussian density distribution (10), also centered on z = 0. The distribution is 
cigar- shaped ellipsoid, with eigenvectors of Q n giving the orientation of various axes at 
time n, see figure 1 1 (b). If the fixed point is attractive, with all multipliers of M strictly 
contracting, any compact initial measure (not only initial distributions of Gaussian form) 
converges under applications of (14) to the unique invariant natural measure po(z) whose 
covariance matrix satisfies the condition 

Q = MQM T + A . (17) 

For a repelling fixed point the condition (15) on the adjoint eigenvector po yields 

MQM T = Q + A, (18) 

with a very different interpretation: as the Jacobian matrix M has only expanding Flo- 
quet multipliers, the deterministic dynamics expands the fixed-point neighborhood ex- 
ponentially, with no good notion of a local neighborhood in the large forward time limit. 
Instead, its past defines the neighborhood, with Q the covariance of the optimal distribu- 
tion of points that can reach the fixed point in one time step, given the diffusion tensor 
A. 

These conditions are central to control theory, where the attracting fixed point con- 
dition (17) is called the Lyapunov equation (see appendix C), Q and Q are known re- 
spectively as controllability and observability Gramians, and there is much wisdom and 
open source code available to solve these (see remark C.l), as well as the more general 
hyperbolic equations. In order to develop some intuition about the types of solutions we 
shall encounter, we assume first, for illustrative purposes, that [dxd] Jacobian matrix M 
has distinct real contracting Floquet multipliers {Ai, A2, • • • , A^} and right eigenvectors 
Me^ = Aje^ . Construct from the d column eigenvectors a [dxd] similarity transfor- 
mation 

S=(e(V 2 ),--.,eW) 

that diagonalizes M, S~ l MS = A and its transpose S T M T (S~ l ) T = A. Define Q = 
S^QiS- 1 ) 7 and A = 5- 1 A(5 _1 ) r . The fixed point condition (17) now takes form 
Q — AQA = A . The matrix elements are Qij{ \ — A ; Aj) = A ;J - , so 



and the attracting fixed point covariance matrix in the original coordinates is given by 



Q = SQS T . (20) 
For the adjoint case, the same algebra yields 



^•=a,.a;-i- (2i) 



for the matrix elements of Q, with the covariance matrix in the fixed coordinates again 
given by Q = SQS T . 

As (20) is not a similarity transformation, evaluation of the covariance matrix Q 
requires a numerical diagonalization, which yields the singular values and singular 
vectors (principal axes) of the equilibrium Gaussian 'cigar' (see appendix A). The 
singular vectors of this symmetric matrix have their own orientations, distinct from the 
left/right eigenvectors of the non-normal Jacobian matrix M. 

Remark 3.2 Hyperbolic flows. The methods to treat the cases where some of the eigen- 
directions are unstable, and other stable are implicit in the Oseledec [30] definition of Lyapunov 
exponents, the rigorous proof of existence of classical spectral (Fredholm) determinants by 
Rugh [31], and the controllability and observability Gramians of control theory [32]: the flow at a 
hyperbolic fixed point or cycle point can be locally factorized into stable and unstable directions, 
and for unstable directions one needs to study noise evolution in the past, by means of the adjoint 
operator (6). 



3.2. In nonlinear world noise is never isotropic 

Now that we have established the exact formulas (14), (15) for the extent of the noise- 
smeared out neighborhood of a fixed point, we turn to the problem of computing them 
for periodic orbits. An attractive feature of the deterministic periodic orbit theory is that 
certain properties of periodic orbits, such as their periods and Floquet multipliers, are in- 
trinsic, independent of where they are measured along the the periodic orbit, and invari- 
ant under all smooth conjugacies, i.e., all smooth nonlinear coordinate transformations. 
Noise, however, is specified in a given coordinate system and breaks such invariances 
(for an exception, a canonically invariant noise, see Kurchan [33]). Each cycle point has 
a different memory and differently distorted neighborhood, so we need to compute the 
Fokker-Planck eigenf unction p a at each cycle point x a . 

The basic idea is simple: A periodic point of an n-cycle is a fixed point of the nth 
iterate of the map (1). Hence the formula (16) for accumulated noise, together the fixed 
point condition (17) also yields the natural measure covariance matrix at a periodic point 
x a on a periodic orbit p, 



Qa = Mp^Q a Mp a + A pja , 



(22) 



where 



A Pja = Aa + Ma-iAa-iM^+M^Aa^iM^f 

+ ■■ +K P - n l p+1 ^-n p+ l (M n a r n l p+1 ) T (23) 

is the noise accumulated per a single transversal of the periodic orbit, M pA = M p {x a ) is 
the cycle Jacobian matrix (70) evaluated on the periodic point x a , and we have used the 
periodic orbit condition x a+np = x a . Similarly, for the adjoint evolution the fixed point 
condition (18) generalizes to 

M p , a Q a Ml a = Qa + Ap,a, (24) 

where 

A p , a = A a + M a+l A a+l Ml +l +Mj +2 A a+2 (Mj +2 ) T 

+--+<CA + n p -i(M;;;;_ 1 ) r (25) 

is the noise accumulated per a single transversal of the periodic orbit backward in time. 

As there is no single coordinate frame in which different M k a _ k A a ^{M k a _ k ) T can be 
simultaneously diagonalized, the accumulated noise is never isotropic. So the lesson is 
that regardless of whether the external noise A is isotropic or anisotropic, the nonlinear 
flow always renders the effective noise anisotropic and spatially inhomogeneous. 



4. ONE-DIMENSIONAL INTUITION 

The very general, exact formulas that we have obtained so far (and so easily), valid in 
any dimension, might be elegant, but it is a bit hard to get one's head around a formula 
such as the expression for the accumulated cycle noise (23). These results are easier to 
grasp by studying the effect of noise on 1 -dimensional systems, such as the noisy linear 
map (7), 

Zn+l = f(z n ) + £n , f(z n ) = Az„ , (26) 

with the deterministic fixed point at f(z) = z = 0, and additive white noise (2) with 
variance A. The density p (x) of trajectories evolves by the action of the Fokker- Planck 
operator (3): 

[J2fp](x) = j[dy]e-l^p{y). (27) 

If a 1 -dimensional noisy linear map (26) is contracting, any initial compact measure 
converges under applications of (27) to the unique invariant natural measure po(z) 
concentrated at the deterministic fixed point z = whose variance (10) is given by (17): 

e=i^' ^=vm e ~ zl/2Q - (28) 

The variance (22) of a periodic point x a on an attractive n -cycle p is 

G «=r=A2' A=/: '' (29) 



where the accumulated noise per a cycle traversal (23) is given by 



A P ,a = a (i + (f^f + {fU? +■■■+ (/r,;; +1 ) 2 ) . oo) 

Variance (28) expresses a balance between contraction by A and diffusive smearing by 
A at each time step. For strongly contracting A, the width is due to the noise only. As 
|A| — > 1 the width diverges: the trajectories are only weakly confined and diffuse by 
Brownian motion into a broad Gaussian. 

Consider next the adjoint operator acting on a repelling noisy fixed point, |A| > 1. 
The stationary measure condition (18) yields 

While the dominant feature of the attracting fixed point variance (28) was the diffusion 
strength A, weakly modified by the contracting multiplier, for the unstable fixed point 
the behavior is dominated by the expanding multiplier A; the more unstable the fixed 
point, the smaller is the neighborhood one step in the past that can reach it. 
The variance (24) of a periodic point x a on an unstable n -cycle p is 

& = T^{W- + (ft?-'-^ ] - (32) 

For an unstable cycle typically all derivatives along the cycle are expanding, \f' a+k \ > 1, 
so the dominant term in (32) is the most recent one, Q a ps A/ (fa) 2 . By contrast, forward 
in time (29) the leading estimate of variance of an attractive periodic point is Q a ~ A. 
These leading estimates are not sensitive to the length of the periodic orbit, so all 
trajectories passing through a neighborhood of periodic point x a will have comparable 
variances. 




4.1. Ornstein-Uhlenbeck spectrum 

The variance (17) is stationary under the action of ££ , and the corresponding Gaus- 
sian is thus an eigenf unction. Indeed, as we shall now show, for the linear flow (27) 
the entire eigenspectrum is available analytically, and as Q a can always be brought to 
a diagonal, factorized form in its orthogonal frame, it suffices to understand the sim- 
plest case, the Ornstein-Uhlenbeck process (see appendix B.l) in one dimension. The 
linearized Fokker-Planck operator is a Gaussian, so it is natural to consider the set of 
Hermite polynomials, Hq(x) = 1, H\(x) = 2x, ^(x) = Ax 2 — 2, • • • , as candidates for 
its eigenfunctions. H n (x) is an nth-degree polynomial, orthogonal with respect to the 
Gaussian kernel 

nn ] n [ dx H m (x) e- x2 H n (x) = 8 mn . (33) 
l n n \ y TZ J 

There are three cases to consider: 



|A| > 1 expanding case: The form of the left po eigenfunction (31) suggests that we 
rescale x — > xj \/2Q and absorb the Gaussian kernel in (33) into left eigenfunctions po, 
Pi,---, 

fcW = ^34^ ff ' ((2r ' /2;> '" W ' <34) 
The right eigenfunctions are then 

p k (z) = (2Q) k / 2 H k ((2Q)- l / 2 z), (35) 
By construction the left, right eigenfunctions are orthonormal to each other: 

j dx p k (x) p j (x) = 8 kj . (36) 

One can verify [3] that for the fixed point z = 0, these are the right, left eigenfunctions 
of the adjoint Fokker- Planck operator (6), where the Mi eigenvalue is l/|A]A fc . Note 
that the Floquet multipliers A* are independent of the noise strength, so they are the 
same as for the A — > deterministic Perron-Frobenius operator (75). 

|A| = 1 marginal case: This is the pure diffusion limit, and the behavior is not 
exponential, but power-law. If the map is nonlinear, one needs to go to the first non- 
vanishing nonlinear order in Taylor expansion (71) to reestablish the control [12]. This 
we do in sect. 7. 

|A| < 1 contracting case: In each iteration the map contracts the cloud of noisy 
trajectories by Floquet multiplier A toward the x = fixed point, while the noise smears 
them out with variance A. Now what was the left eigenfunction for the expanding case 
(34) is the peaked right eigenfunction of the Fokker-Planck operator, {po, pi, pi,- ■ ■ }, 
with eigenvalues {1, A, A 2 ,- ■ • } [11, 12] 

p k (x) = N k l H k {{2Q)- l l 2 x)e-* 2 l 2 Q, £ = A/(1-A 2 ), (37) 

where H k (x) is the Mi Hermite polynomial, and N k l follows from the prefactor in (34). 

These discrete time results can be straightforwardly generalized to continuous time 
flows of sect. B, as well as to higher dimensions. So far we have used only the leading 
eigenfunctions (the natural measure), but in sect. 6 we shall see that knowing the whole 
spectrum in terms of Hermite polynomial is a powerful tool for the computation of 
weak-noise corrections. 

Remark 4.1 Ornstein-Uhlenbeck process. The simplest example of a continuous time 
stochastic flow (86) is the Langevin flow (99) in one dimension. In this case, nothing is lost by 
considering discrete-time dynamics which is strictly equivalent to the continuous time Ornstein- 
Uhlenbeck process (100) discussed in appendix B. 1. 
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FIGURE 3. (a) /o,/i: branches of the deterministic map (38) for Ao = 8 and b = 0.6. The local 
eigenfunctions p a o with variances given by (32) provide a state space partitioning by neighborhoods of 
periodic points of period 3. These are computed for noise variance A = 0.002. The neighborhoods ^#ooo 
and ^#ooi already overlap, so cannot be resolved further, (b) The next generation of eigenfunctions 
shows how the neighborhoods of the optimal partition cannot be resolved further. For periodic points of 
period 4, only ^#011 can be resolved further, into ./#ono an d -^bm (second and third peak from the left), 
but that would not change the transition graph of figure 5. 



5. HAVE NO FEAR OF GLOBALIZATION 

We are now finally in position to address our challenge: Determine the finest possible 
partition for a given noise. 

We shall explain our 'the best possible of all partitions ' hypothesis by formulating 
it as an algorithm. For every unstable periodic point x a of a chaotic one-dimensional 
map, we calculate the corresponding width Q a of the leading Gaussian eigenfunction of 
the local adjoint Fokker- Planck operator Jzf^. Every periodic point is assigned a one- 
standard deviation neighborhood [x a — \/Qa,x a + \Z~Qa]- We cover the state space with 
neighborhoods of orbit points of higher and higher period n p , and stop refining the local 
resolution whenever the adjacent neighborhoods, say of x a and x^, overlap in such a 
way that \x a — X\,\ < v2a + \/~Qb- As an illustration of the method, consider the chaotic 
repeller on the unit interval 

x n+ i =Aox„(l-x n )(l-bx n ) + £ n , A = 8, 6 = 0.6, (38) 

with noise strength A = 0.002. 

The map is plotted in figure 3 (a), together with the local eigenfunctions p a with vari- 
ances given by (32). Each Gaussian is labeled by the {/o ? /i} branches visitation se- 
quence of the corresponding deterministic periodic point (a symbolic dynamics, how- 
ever, is not a prerequisite for implementing the method). Figure 3 (b) illustrates the 
overlapping of partition intervals: {^#0007*^001 }> {^0101,^0100} overlap and so do 
all other neighborhoods of the period n p = 4 cycle points, except for ^#0110 an d ^blll- 
We find that in this case the state space (the unit interval) can be resolved into 7 neigh- 
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FIGURE 4. (upper panel) The unit interval partitioned deterministic ally by a binary tree. Due to the 
noise, the partitioning stops where the eigenfunctions of figure 3 overlap significantly, (lower panel) Once 
the optimal partition is found, the symbolic dynamics is recoded by relabeling the finite partition intervals, 
and refashioned into the transition graphs of figure 5. 



borhoods 

{^#oo, -^on, ^oio, ^no, «#m, ^101, ^ioo}- (39) 

It turns out that resolving ^#011 further into ^ono and ^oin would not affect our 
estimates, as it would produce the same transition graph. 

Once the finest possible partition is determined, a finite binary tree like the one 
in figure 4 is drawn: Evolution in time maps the optimal partition interval — * 
{^#110, ^in}, ^#oo — > {^00)^011)^010}* etc.. This is summarized in the transition 
graph in figure 5, which we will use to estimate the escape rate and the Lyapunov 
exponent of the repeller. 

Remark 5.1 A brief history of state space partitions. There is considerable prior literature 
that addresses various aspects of the 'optimal partition' problem. Before reviewing it, let us state 
what is novel about the optimal partition hypothesis formulated here: Our estimates of limiting 
resolution are local, differing from region to region, while all of the earlier limiting resolution 
estimates known to us are global, based on global averages such as Shannon entropy or quantum- 
mechanical H 'granularity' of phase space. We know of no published algorithm that sets a limit 
to the resolution of a chaotic state space by studying the interplay of the noise with the local 
stretching/contracting directions of the deterministic dynamics, as we do here. 

The engineering literature on optimal experimental design [34, 35, 36, 37] employs criteria 
such as 'D-optimality,' the maximization of the Shannon information content of parameter esti- 
mates. Purely statistical in nature, these methods have little bearing on the dynamical approach 
that we pursue here. 

In 1983 Crutchfield and Packard [38] were the first to study the problem of an optimal 
partition for a chaotic system in the presence of noise, and formulate a state space resolution 
criterion in terms of a globally averaged "attainable information." The setting is the same that we 
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FIGURE 5. (a) Transition graph (graph whose links correspond to the nonzero elements of a transition 
matrix Tt, a ) describes which regions b can be reached from the region a in one time step. The 7 nodes 
correspond to the 7 regions of the optimal partition (39). Dotted links correspond to symbol 0, and the full 
ones to 1, indicating that the next region is reached by the /o, respectively f\ branch of the map plotted 
in figure 3. (b) The region labels in the nodes can be omitted, with links keeping track of the symbolic 
dynamics. 



assume here: the laws governing deterministic dynamics are given, and one studies the effects 
of noise (be it intrinsic, observational or numerical) on the dynamics. They define the most 
efficient symbolic encoding of the dynamics as the sequence of symbols that maximizes the 
metric entropy of the entire system, thus their resolution criterion is based on a global average. 
Once the maximum for a given number of symbols is found, they refine the partition until the 
entropy converges to some value. They formulate their resolution criterion in terms of attainable 
information, a limiting value for the probability to produce a certain sequence of symbols from 
the ensemble of all possible initial conditions. Once such limit is reached, no further refinements 
are possible. 

Most of the dynamical systems literature deals with estimating partitions from observed 
data [39]. Tang and co-workers [40] assume a noisy chaotic data set, but with the laws of 
dynamics assumed unknown. Their method is based on maximizing Shannon entropy and at 
the same time minimizing an error function with respect to the partition chosen. The same 
idea is used by Lehrman et al. [41] to encode chaotic signals in higher dimensions, where they 
also detect correlations between different signals by computing their conditional entropy. For a 
review of symbolic analysis of experimental data up to 2001, see Daw, Finney and Tracy [39]. 

Kennel and Buhl [42, 43, 44] estimate partitions for (high-dimensional) flows from noisy 
time-series data by minimizing a cost function which maximizes the correlation between dis- 
tances in the state space and in the symbolic space, and indicates when to stop adjusting their 
partitions and therefore what the optimal partition is. In ref. [42] their guiding principle for 
a good partition is that short sequences of consecutive symbols ought to localize the corre- 
sponding continuous state space point as well as possible. They embed symbol sequences into 
the unit square, and minimize the errors in localizing the corresponding state space points un- 
der candidate partitions. Holstein and Kantz [45] present an information-theoretic approach to 
determination of optimal Markov approximations from time series data based on balancing the 
modeling and the statistical errors in low-dimensional embedding spaces. Boland, Galla and 
McKane [46] study the effects of intrinsic noise on a class of chemical reaction systems which 
in the deterministic limit approach a limit cycle in an oscillatory manner. 



A related approach to the problem of the optimal resolution is that of the refinement of a 
transition matrix: given a chaotic, discrete-time dynamical system, the state space is partitioned, 
and the probabilities of points mapping between regions are estimated, so as to obtain a transition 
matrix, whose eigenvalues and eigenfunctions are then used to evaluate averages of observables 
defined on the chaotic set. The approach was first proposed in 1960 by Ulam [47, 2], for 
deterministic dynamical systems. He used a uniform-mesh grid as partition, and conjectured that 
successive refinements of such coarse-grainings would provide a convergent sequence of finite- 
state Markov approximations to the Perron-Frobenius operator. Rechester and White [48, 49] 
have proposed dynamics-based refinement strategies for constructing partitions for chaotic maps 
in one and two dimensions that would improve convergence of Ulam's method. 

Bollt et al. [50] subject a dynamical system to a small additive noise, define a finite Markov 
partition, and show that the Perron-Frobenius operator associated to the noisy system is repre- 
sented by a finite-dimensional stochastic transition matrix. Their focus, however, is on approx- 
imating the natural measure of a deterministic dynamical system by the vanishing noise limit of 
a sequence of invariant measures of the noisy system. 

In ref. [44] Kennel and Buhl approximate the distribution of the points in each symbolic 
region by a Gaussian with mean jx and variance t, 



/(x| ^' T)= (2^ eXP 



(40) 



and estimate "code length" by the ad hoc Rissanen prior on x, defined with no reference to 
dynamics, and thus morally unrelated to our periodic orbits based optimal partition. 

Dellnitz and Junge [53], Guder and Kreuzer [54], Froyland [55], and Keane et al. [56] propose 
a variety of non-uniform refinement algorithms for such grids, reviewed in a monograph by 
Froyland [57], who also treats their extension to random dynamical systems. In all cases, the 
ultimate threshold for every refinement is determined by the convergence of the spectrum of the 
transition matrix. 

Theoretical investigations mostly focus on deterministic limits of stochastic models. The 
Sinai-Ruelle-Bowen [65, 66, 67] or natural measure (also called equilibrium measure, SRB 
measure, physical measure, invariant density, natural density, or natural invariant) is singled 
out amongst all invariant measures by its robustness to weak-noise perturbations, so there is 
considerable literature that studies it as the deterministic limit of a stochastic process. 



6. FINITE FOKKER-PLANCK OPERATOR, AND STOCHASTIC 

CORRECTIONS 

Next we show that the optimal partition enables us to replace Fokker-Planck PDEs by 

in 

finite-dimensional matrices. The variance (32) is stationary under the action of ££ a p , 
and the corresponding Gaussian is thus an eigenfunction. Indeed, as we showed in 
sect. 4.1, for the linearized flow the entire eigenspectrum is available analytically. For 
a periodic point x a E p, the n p th iterate J?f a p of the linearization (9) is the discrete time 
version of the Ornstein-Uhlenbeck process, with left po, pi, • • • , respectively right po, 
pi, • • ■ mutually orthogonal eigenfunctions (34). 

The number of resolved periodic points determines the dimensionality of the Fokker- 
Planck matrix. Partition (39) being the finest possible partition, the Fokker-Planck oper- 



ator now acts as [7 x 7] matrix with non-zero a — >■ b entries expanded in the Hermite 
basis, 

[Ua]kj = (pb,k\^\Paj) 

= r dz b dzaP c _ ( is Zh) 2_ (z b -f^a)) 2 
J 2j+ l j\n^X/l 

xH k (j5z b )Hjil5za), (41) 

where 1//3 = ^/2Qa, and z a is the deviation from the periodic point x a . 

Periodic orbit theory (summarized in appendix A) expresses the long-time dynamical 
averages, such as Lyapunov exponents, escape rates, and correlations, in terms of the 
leading eigenvalues of the Fokker-Planck operator. In our optimal partition approach, Jzf 
is approximated by the finite-dimensional matrix L, and its eigenvalues are determined 
from the zeros of det(l — zL), expanded as a polynomial in z, with coefficients given by 
traces of powers of L. As the trace of the nth iterate of the Fokker-Planck operator ££ n 
is concentrated on periodic points f n (x a ) = x a ,we evaluate the contribution of periodic 
orbit p to trL"^ by centering L on the periodic orbit, 

t p = tr p = trL ad ■ ■ ■ L ch L ba , (42) 

where x a ,x b ,---x c i E p are successive periodic points. To leading order in the noise 
variance A, t p takes the deterministic value t p = l/|A p — 1|. The nonlinear diffusive 
effects in (41) can be accounted for [19] by the weak- noise Taylor series expansion 
around the periodic point x a , 

e =— =e =— (l-V2A(f a f a zl + f a z 2 a z b ) + 0(A)). (43) 

Such higher order corrections will be needed in what follows for a sufficiently accurate 
comparison of different methods. 

We illustrate the method by calculating the escape rate y = — lnzo, where z 1 is the 
leading eigenvalue of Fokker-Planck operator Jzf , for the repeller plotted in figure 3. 
The spectral determinant can be read off the transition graph of figure 5 and its loop 
expansion in figure 6, 

det( 1 - zL) = 1 - (to + h)z - (/01 - toh ) z 2 

— (tooi +ton —toito — toiti)z 3 

— (?ooi 1 +*om — fooi^i — *on*o — *on*i + *oi*o*i)z 
-(^ooin -^oiii^o-^ooii^i +tontoti)z 5 

— (^ooion + foonoi — foon^oi — fooi?on)z 6 

— (^0010111 +^0011101 — *ooion*i — *oonoi*i 

— *oom*oi + *oon*oi*i + ^001^011^1) z 7 - (44) 

The polynomial coefficients are given by products of non-intersecting loops of the tran- 
sition graph [24], with the escape rate given by the leading root z 1 of the polynomial. 
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FIGURE 6. (a)-(i) The fundamental cycles for the transition graph figure 5 (b), i.e., the set of its non- 
self-intersecting loops. Each loop represents a local trace t p ; together they form the determinant (44). 



Twelve p eriodic orbits 0, 1,01,001,011,0011,0111,00111,001101,001011,0010111, 
0011101 up to period 7 (out of the 41 contributing to the noiseless, deterministic cycle 
expansion up to cycle period 7) suffice to fully determine the spectral determinant of the 
Fokker-Planck operator. In the evaluation of traces (42) we include stochastic correc- 
tions up to order O(A) (an order beyond the term kept in (43)). The escape rate of the 
repeller of figure 3 so computed is reported in figure 7. 

Since our optimal partition algorithm is based on a sharp overlap criterion, small 
changes in noise strength A can lead to transition graphs of different topologies, and it is 
not clear how to assess the accuracy of our finite Fokker-Planck matrix approximations. 
We make three different attempts, and compute the escape rate for: (a) an under-resolved 
partition, (b) several deterministic, over-resolved partitions, and (c) a direct numerical 
discretization of the Fokker-Planck operator. 

(a) In the example at hand, the partition in terms of periodic points 00, 01, 11 and 10 
is under-resolved; the corresponding escape rate is plotted in figure 7. (b) We calculate 
the escape rate by over-resolved periodic orbit expansions, in terms of all deterministic 
periodic orbits of the map up to a given period, with t p evaluated in terms of Fokker- 
Planck local traces (42), including stochastic corrections up to order O(A). Figure 7 
shows how the escape rate varies as we include all periodic orbits up to periods 2 through 
8. Successive estimates of the escape rate appear to converge to a value different from 
the optimal partition estimate, (c) Finally, we discretize the Fokker-Planck operator 
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FIGURE 7. (a) The escape rate y of the repeller in figure 3 plotted as function of number of partition 
intervals TV, estimated using: (♦) under-resolved 4-interval and the 7-interval optimal partition, (•) all 
periodic orbits of periods up to n = 8 in the deterministic, binary symbolic dynamics, with TV, = 2" 
periodic -point intervals (the deterministic, noiseless escape rate is y<> = 0.7011), and (■) a uniform 
discretization (45) in/V = 16, • • • ,256 intervals. For /V = 512 discretization yields y num = 0.73335(4). (b) 
Number of neighborhoods required by the optimal partition method vs. the noise strength A. 



by a piecewise-constant approximation on a uniform mesh on the unit interval [47], 

[Sf]tj = n^^= / dx ! dye-^-m)\ (45) 

where ^ is the z'th interval in equipartition of the unit interval into /V equal segments. 
Empirically, N = 128 intervals suffice to compute the leading eigenvalue of the dis- 
cretized [128 x 128] matrix [J5f]y to four significant digits. This escape rate, figure 7, is 
consistent with the N = 1 optimal partition estimate to three significant digits. 

We estimate the escape rate of the repeller (38) for a range of values of the noise 
strength A. The optimal partition method requires a different numbers of neighborhoods 
every time for different noise strengths. The results are illustrated by figure 7 (b) and 
8, with the estimates of the optimal partition method within 2% of those given by the 
uniform discretization of Fokker-Planck. One can also see from the same table that the 
escape rates calculated with and without higher order corrections to the matrix elements 
(41) are consistent within less than 2%, meaning that the stochastic corrections (43) do 
not make a significant difference, compared to the effect of the optimal choice of the 
partition, and need not be taken into account in this example. 

The optimal partition estimate of the Lyapunov exponent is given by A = 
(In |A|) / (n), where the cycle expansion average of an integrated observable A [24] 

(A) = A t +A 1 ti+ [Aoi?oi-(A +Ai)?u*i] 
+ [Aooi^ooi - (A i +A )t it ] + ■■■ 
+ [Aoii*oii - (Aoi +Ai>oi*i] + • • • (46) 

is the finite sum over cycles contributing to (44), and In \A p \ = £ln |/'(;c a )|, the sum over 
the points of cycle p, is the cycle Lyapunov exponent. On the other hand, we also use 
the discretization (45) to cross check our estimate: this way the Lyapunov exponent is 
evaluated as the average 

A = J dxe r p(x)ln\f(x)\, (47) 
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FIGURE 8. (a) Escape rates of the repeller (38) vs. the noise strength A, using: the optimal partition 
method with (♦) and without (x) stochastic corrections; (■) a uniform discretization (45) in N = 128 
intervals, (b) The Lyapunov exponent of the repeller (38) vs. the noise strength A, using: the optimal 
partition method (•) without stochastic corrections, and (o) a uniform discretization (45) over AT = 128 
intervals. 



where p(x) is the leading eigenfunction of (45), y is the escape rate, and e r p is the 
normalized repeller measure, / dxe r p(x) = 1. Figure 8 shows close agreement (< 1%) 
between the Lyapunov exponent estimated using the average (46), where t p = l/\A p — 1 1 
(no higher-order stochastic corrections), and the same quantity evaluated with (47), by 
the discretization method (45). 

Remark 6.1 Weak noise corrections. The weak-noise corrections to the spectrum of 
evolution operators were first treated by Gaspard [4] for continuous-time systems, and in a 
triptych of articles [19, 20, 21] for discrete-time maps: they can be computed perturbatively 
to a remarkably high order [69] in the noise strength A. However, as we have shown here, the 
eigenvalues of such operators offer no guidance to the 'optimal partition' problem; one needs to 
compute the eigenf unctions. 



7. WHEN THE GAUSSIAN APPROXIMATION FAILS 

The state space of a generic deterministic flow is an infinitely interwoven hierarchy of 
attracting, hyperbolic and marginal regions, with highly singular invariant measures. 
Noise has two types of effects. First, it feeds trajectories into state space regions that 
are deterministically either disconnected or transient ("noise induced escape," "noise 
induced chaos") and second, it smoothens out the natural measure. Here, we are mostly 
concerned with the latter. Intuitively, the noisy dynamics erases any structures finer 
than the optimal partition, thus -in principle- curing both the affliction of long-period 
attractors/elliptic islands with very small immediate basins of attraction/ellipticity, and 
the slow, power-law correlation decays induced by marginally stable regions of state 
space. So how does noise regularize nonhyperbolic dynamics? 

As a relatively simple example, consider the skew Ulam map [70], i.e., the cubic map 
(38) with the parameter Ao = 1 / f{x c ) ■ The critical point x c is the maximum of / on the 
unit interval, with vanishing derivative f'(x c ) = 0. As this map sends the unit interval 
into itself, there is no escape, but due to the quadratic maximum the (deterministic) 



natural measure exhibits a spike (xh — x)~ l l 2 near the critical value f(x c ) = x^ (see, 
for example, ref. [71] for a discussion). As explained in ref. [70], a close passage to 
the critical point effectively replaces the accumulated Floquet multiplier by its square 
root. For example, for the skew Ulam map (38) n-cycles whose itineraries are of form 
0" _1 1 spend long time in the neighborhood of xq = 0, and then pass close to x c . In the 
neighborhood of xq = the Floquet multiplier gains a factor ~ Ao = f'(xo) for each of 
the first n— 1 iterations, and then experiences a strong, square root contraction during the 

nil 

close passage to the critical point x c , resulting in the Floquet multiplier Ao-oi 06 A 7 , 
and a Lyapunov exponent that converges to Ao/2, rather than Ao that would be expected 
in a hyperbolic flow for a close passage to a fixed point xq. The same strong contraction 
is experienced by the noise accumulated along the trajectory prior to the passage by 
the critical point, rendering, for example, the period-doubling sequences more robust to 
noise than one would naively expect [72, 73, 16]. 

For the corresponding noisy map (1) the critical point is extended into the 'flat top' 
region where \ f'(x) \ <C 1, and the linearized, Gaussian approximation (9) to the Fokker- 
Planck operator does not hold. Thus, we should first modify our choice of densities and 
neighborhoods, as the whole construction leading to the optimal partition algorithm was 
based on the Gaussian approximation. 

The adjoint Fokker-Planck operator acts on a Gaussian density centered at x a , as in 
(10): 

[^p a ](x) = if^e-M^^] 

J — 00 
1 (f(x)-xa) 2 

= e e«+ A . (48) 

C a -1 

Suppose the point x a -\ = f~ l (x a ) around which we want to approximate the new 
density, is very close to the critical point, so that we can write 

p^x) = —e~ 2 (e«+A) = —e'fa-i *S-i/8(G«+A) . ( 49 ) 
C a -\ C a -\ 

During a close passage to the critical point, the variance does not transform linearly, but 
as a square root: 

n _ AV /: -' ¥/8(ga+A W _ W) / 8( 6a + A) \ V2 
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We show in appendix D that in the next iteration the variance of the density p a _i (z a -i) 
transforms again like the variance of a Gaussian, up to order 0(A) in the noise strength. 
By the same procedure, one can again assume the next preimage of the map x a -3 is such 
that the linear approximation is valid, and transform the density p a -i{za-2) (Eq. (115), 
appendix D) up to 0(A) and obtain the same result for the variance, that is 

n Qa-l+^+f'll) ( , u 

Qa-3 = ^2~?2^ ( 51 ) 

Ja-2Ja-3 



which is again the evolution of the variances in the Gaussian approximation. In other 
words, the evolution of the variances goes back to be linear, to 0(A), although the 
densities transformed from the 'quartic Gaussian' (49) are no longer Gaussians. 

The question is now how to modify the definition of neighborhoods given in sect. 5, 
in order to fit the new approximation. Looking for eigenfunctions of j£f' seems to be a 
rather difficult task to fulfill, given the functional forms (49) and (115) involved. Since 
we only care about the variances, we define instead the following map 



Qa-l 




i/: 2 -!<ii 

otherwise 



(52) 



C = 2v / 2r(3/4)/r(l/4), for the evolution of the densities, and take its periodic points 
as our new neighborhoods. In practice, one can compute these numerically, but we will 
not need orbits longer than length n p = 4 in our tests of the partition, therefore we can 
safely assume only one periodic point of f(x) to be close to the flat top, and obtain 
analytic expressions for the periodic points of (52): 
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(53) 



7 



with A p = fa-n+ifa-i 1 ' * s vau d when the cycle starts and ends at a point x a close to the 
flat top. Otherwise, take the periodic point x a ^, that is the k— th preimage of the point 
x a . The corresponding periodic point variance has the form 



Qa-k * tJ-^ (A(l + ... + {ftlf) + Qa) 



,) 2 



(54) 



both expressions (53) and (54) are approximate, as we further assumed AA 2 3> 1, which 

is reasonable when A e [10~ 4 , 10~ 2 ], our range of investigation for the noise strength. 

As before, a neighborhood of width [x a — \Z~Q a ,x a + \Z3 a ] i s assigned to each periodic 
point x a , and an optimal partition follows. However, due to the geometry of the map, 
such partitions as 



{.#000i [^#001,^011] ,-^010,^110,^111,^10} 



(55) 



can occur. In this example the regions -#001 and ^#011 overlap, and the partition results 
in a transition graph with three loops (cycles) of length one, while we know that our map 
only admits two fixed points. In this case we decide to follow the deterministic symbolic 
dynamics and ignore the overlap. 
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FIGURE 9. (a) Escape rate y of the 'skew Ulam' map vs. noise strength A, using: (x) the optimal 
partition method; (■) a uniform discretization (45) in A = 128 intervals, (b) Number of neighborhoods 
required by the optimal partition method vs. the noise strength A. 



Let us now test the method by estimating once again the escape rate of the noisy map 
(38). We note that the matrix elements 

[Lfajjfcy = (pb,k\-&\Pa,j) 

_ f dZbdZaP _^ Zb )2_(z b -^a) 2 

J Vj\7tV2A e 

xH k (j5z b )Hj(l5za), (56) 

1/j8 = \J1Q a , should be redefined in the neighborhood of the critical point of the map, 
where the Gaussian approximation to Jz? fails. We follow the approximation made in 
(49): 

\ ^^-"^^m'iSB^ (57) 
J 2 J ji7tV2A 

However, as A decreases, it also reduces the quadratic term in the expansion of the 
exponential, so that the linear term f' a z a must now be included in the matrix element: 

f dZhdZaP ffl, ^ ^b-f'aza-f'a^lm 2 ,„ . 

[Ua]kj = / * %r e- (Pzb) s H k ((5z b )H0z a ) , (58) 

J 2 J ji7ZV2A 

We find in our model that the periodic orbits we use in our expansion have jc a 's within 
the flat top, such that f a ~ 10 _1 and f a ~ 10, and therefore (57) better be replaced 
with (58) when A ~ 10 -4 . In order to know whether a cycle point is close enough to 
the flat top for the Gaussian approximation to fail, we recall that the matrix element 
(56) is the zeroth-order term of a series in A, whose convergence can be probed by 
evaluating the higher order corrections (43): when the 0(y/A) and O(A) corrections 
are of an order of magnitude comparable or bigger than the one of (56), we conclude 
that the Gaussian approximation fails and we use (57) or (58) instead. Everywhere else 
we use our usual matrix elements (56), without the higher-order corrections, as they 
are significantly larger than in the case of the repeller, and they are not accounted for 



by the optimal partition method, which is entirely based on a zeroth-order Gaussian 
approximation of the evolution operator. Like before, we tweak the noise strength A 
within the range [10~ 4 , 10~ 2 ] and compare the escape rate evaluated with the optimal 
partition method and with the uniform discretization (45). The results are illustrated in 
figure 9: the uniform discretization and the method of the optimal partition are consistent 
within a 5% margin. 

8. SUMMARY AND CONCLUSIONS 

Physicists tend to believe that with time Brownian motion leads to x(t) 2 ^At broadening 
of a noisy trajectory neighborhood. In nonlinear dynamics nothing of the sort happens; 
the noise broadening is balanced by non-linear stretching and contraction, and infinite 
length recurrent Langevin trajectories have finite noise widths, not widths that spread oc 
time. Here periodic orbits play special role: computable in finite time, they persist for for 
infinite time, and are thus natural objects to organize state space partitions around. On 
the other hand, computation of unstable periodic orbits in high-dimensional state spaces, 
such as Navier-Stokes, is at the border of what is currently feasible numerically [74, 75], 
and criteria to identify finite sets of the most important solutions are very much needed. 
Where are we to stop calculating orbits of a given hyperbolic flow? Intuitively, as we 
look at longer and longer periodic orbits, their deterministic neighborhoods shrink ex- 
ponentially with time, while the variance of the noise-induced orbit smearing remains 
bounded; there has to be a turnover time, a time at which the noise-induced width over- 
whelms the exponentially shrinking deterministic dynamics, so that no better resolution 
is possible. Given a specified noise, we need to find, periodic orbit by periodic orbit, 
whether a further sub-partitioning is possible. 

We have described here the optimal partition hypothesis, a method for partitioning the 
state space of a chaotic repeller in presence of weak Gaussian noise first introduced in 
ref. [7]. The key idea is that the width of the linearized adjoint Fokker- Planck operator 
eigenfunction computed on an unstable periodic point x a provides the scale beyond 
which no further local refinement of state space is feasible. This computation enables 
us to systematically determine the optimal partition, the finest state space resolution 
attainable for a given chaotic dynamical system and a given noise. Once the optimal 
partition is determined, we use the associated transition graph to describe the stochastic 
dynamics by a finite dimensional Fokker- Planck matrix. An expansion of the Fokker- 
Planck operator about periodic points was already introduced in refs. [19, 20, 21], 
with the stochastic trace formulas and determinants [19, 4] expressed as finite sums, 
truncated at orbit periods corresponding to the local turnover times. A novel aspect of 
the work presented here is its representation in terms of the Hermite basis (sect. B.l), 
eigenfunctions of the linearized Fokker-Planck operator (9), and the finite dimensional 
matrix representation of the Fokker-Planck operator. 

It should be noted that our linearization of Fokker-Planck operators does not im- 
ply that the nonlinear dynamics is being modeled by a linear one. Our description is 
fully nonlinear, with periodic orbits providing the nonlinear backbone of chaotic dy- 
namics, dressed up stochastically by Fokker-Planck operators local to each cycle. This 
is a stochastic cousin of Gutz wilier 's WKB approximation based semi-classical quanti- 



zation [76] of classically chaotic systems, where in a parallel effort to utilize quantum- 
mechanical h 'graininess' of the quantum phase space to terminate periodic orbit sums, 
Berry and Keating [77] have proposed inclusion of cycles of periods up to a single 
'Heisenberg time'. In light of the stochastic dynamics insights gained here, this pro- 
posal merits a reexamination - each neighborhood is likely to have its own Heisenberg 
time. 

The work of Abarbanel et al. [78, 79, 28, 80] suggest one type of important applica- 
tion beyond the low-dimensional Fokker-Planck calculations undertaken here. In data 
assimilation in weather prediction the convolution of noise variance and trajectory vari- 
ance (14) is a step in the Kalman filter procedure. One could combine the state space 
charts of turbulent flows of ref. [81] (computed in the full 3-dimensional Navier-Stokes) 
with partial information obtained in experiments (typically a full 3-dimensional velocity 
field, fully resolved in time, but measured only on a 2-dimensional disk section across 
the pipe). The challenge is to match this measurement of the turbulent flow with a state 
space point in a 10 5 -dimensional ODE representation, and then track the experimen- 
tal observation to improve our theoretical prediction for the trajectory in the time ahead. 
That would be the absolutely best 'weather prediction' attainable for a turbulent pipe 
flow, limited by a combination of Lyapunov time and observational noise. In our par- 
lance, the 'optimal partition of state space.' 

We have tested our optimal partition hypothesis by applying it to evaluation of the 
escape rates and the Lyapunov exponents of a Id repeller in presence of additive noise. 
In the Id setting numerical tests indicate that the 'optimal partition' method can be as 
accurate as the much finer grained direct numerical Fokker-Planck operator calculations. 
In higher dimensions (and especially in the extreme high dimensions required by fluid 
dynamics stimulations) such direct Fokker-Planck PDE integrations are not feasible, 
while the method proposed here is currently the only implementable approach. 

The success of the optimal partition hypothesis in a one-dimensional setting is encour- 
aging, and use of noise as a smoothing device that eliminates singularities and patholo- 
gies from clusterings of orbits is promising. However, higher-dimensional hyperbolic 
maps and flows, for which an effective optimal partition algorithm would be very use- 
ful, present new, as yet unexplored challenges of disentangling the subtle interactions 
between expanding, marginal and contracting directions; the method has not yet been 
tested in a high-dimensional hyperbolic setting. A limiting factor to applications of the 
periodic orbit theory to high-dimensional problems ranging from fluid flows to chemi- 
cal reactions might be the lack of a good understanding of periodic orbits in more than 
three dimensions, of their stability properties, their organization and their impact on the 
dynamics. 

In summary: Each periodic point owns a cigar, which for a high-dimensional dissipa- 
tive flow is shaped along a handful of expanding and least-contracting directions. The 
remaining large (even infinite!) number of the strongly contracting directions is limited 
by the noise; the cigar always has the dimensionality of the full state space. Taken to- 
gether, the set of overlapping cigars, or the optimal partition, weaves the carpet (of the 
full dimensionality of state space) which envelops the entire 'inertial manifold' explored 
by turbulent dynamics. 
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A. PERIODIC ORBIT THEORY, DETERMINISTIC DYNAMICS 

We offer here a brief review of deterministic dynamics and periodic orbit theory. All 
of this is standard, but needed to set the notation used above. The reader might want to 
consult ref. [24] for further details. 

Though the main applications we have in mind are to continuous flows, for purposes at 
hand it will suffice to consider discrete time dynamics obtained, for example, by reduc- 
ing a continuous flow to mappings between successive Poincare sections, as in figure 10. 
Consider dynamics induced by iterations of a J-dimensional map / : ^# — > ^# , where 
^ C W* is the state space (or 'phase space') of the system under consideration. The 
discrete 'time' is then an integer, the number of applications of a map. We denote the 
Mi iterate of map / by composition 

/*(*) =/(/*"*(*)) , A*)=*- (59) 
The trajectory of x = xq is the finite set of points Xj = P(x), 

{x ,x h x 2 ,...,x k } = j*,/(*),/ 2 (*),..., /*(*)} , (60) 

traversed in time k, and the orbit of x is the subset M- x of all points of ^ that can be 
reached by iterations of /. Here x^ is a point in the J-dimensional state space j& , and the 
subscript k indicates time. While a trajectory depends on the initial point x, an orbit is 
a set invariant under dynamics. The transformation of an infinitesimal neighborhood 
of an orbit point x under the iteration of a map follows from Taylor expanding the 
iterated mapping at finite time k. The linearized neighborhood is transported by the 
[d x d] Jacobian matrix 



X=Xq 



(61) 



(J(x) for Jacobian, or derivative notation M(x) — > Df(x) is frequently employed in the 
literature.) The formula for the linearization of Mi iterate 

M k (xo)=M(x k _ l )---M(x 1 )M(x ), M ij = df i /dx j , (62) 




Source: ChaosBook.org 



FIGURE 10. (a) A Poincare hypersurface 8?, defined by a condition U (x) = 0, is intersected by the 
x(t) orbit at times ti,t2,t3,U, and closes a cycle (xi,X2,Xi,X4), Xk = x(tk) € 8?, of topological length 4 
with respect to this section. The crossing z does not count, as it in the wrong direction, (b) The same orbit 
reduced to a Poincare return map that maps points in the Poincare section 3? as x n +\ = f(x n ) . In this 
example the orbit of x\ is periodic and consists of the four periodic points (x\,X2,xj,X4). 



in terms of unit time steps M follows from the chain rule for functional composition, 

d 



M 2 (x ) = M(xi)M(xq). 



-j fk(xo) 
y=f(*o) dXi 



We denote by the ith. eigenvalue or multiplier of the Jacobian matrix M (xo), and by 
the £th Floquet or characteristic exponent, with real part jU^ and phase co^>: 

A £ = e k ^=e k ^ )+i ^l (63) 

Jacobian matrix M k (xo) and its eigenvalues (Floquet multipliers) depend on the initial 
point xo and the elapsed time k. For notational brevity we tend to omit this dependence, 
but in general 

A = A( = Ai{xQ,k) , X = X^(xo,k), (0 = G)^\xo,k) , • • • etc., 

depend on the traversed trajectory. 

A periodic point (cycle point) x^ belonging to a periodic orbit (cycle) of period n is a 
real solution of 

f n (x k )=f(f(...f(x k )...))=x k , fc = 0,l,2,...,n-l. (64) 

For example, the orbitofxi in figure lOisthe 4-cycle (x i,X2, X3, X4 ). The time-dependent 
n-periodic vector fields, such as the flow linearized around a periodic orbit, are described 
by Floquet theory. Hence we shall refer to the Jacobian matrix M p (x) = M" (x) evaluated 
on a periodic orbit p as the monodromy or Floquet matrix, and to its eigenvalues {A P) i, 
A Pi 2, • • • , A p d} as Floquet multipliers. They are flow-invariant, independent of the 
choice of coordinates and the initial point in the cycle p, so we label them by their p 



label. We number the eigenvalues in order of decreasing magnitude |Ai | > | A.2 j > • • • > 
\A c i \ , sort them into sets {e,m, c} 

expanding: {A} e = {A pJ : \A pJ \ > 1} 

marginal: {A} m = {A pJ : |A pJ | = 1} (65) 
contracting: {A} c = {A p j : \A p j\ < 1}, 

and denote by A p (no jth eigenvalue index) the product of expanding Floquet multipliers 

A p = Y[A p , e . (66) 

e 

The stretching/contraction rates per unit time are given by the real parts of Floquet 
exponents 

/4° = — ln|A P)i |. (67) 
n p 

They can be loosely interpreted as Lyapunov exponents evaluated on the prime cycle p. 
A periodic orbit p is stable if real parts of all of its Floquet exponents are strictly 

negative, $ < 0. If all Floquet exponents are strictly positive, jU^ > ju W(n > 0, the 
periodic orbit is repelling, and unstable to any perturbation. If some are strictly positive, 
and rest strictly negative, the periodic orbit is said to be hyperbolic or a saddle, and 
unstable to perturbations outside its stable manifold. Repelling and hyperbolic periodic 
orbits are unstable to generic perturbations, and thus said to be unstable. If all jl^ = 0, 
the orbit is said to be elliptic, and if ji^ = for a subset of exponents, the orbit is said to 
be partially hyperbolic. If all Floquet exponents (other than the vanishing longitudinal 
exponent) of all periodic orbits of a flow are strictly bounded away from zero, the flow 
is said to be hyperbolic. Otherwise the flow is said to be nonhyperbolic. 

The Jacobian matrix M is in general non-normal: it neither symmetric, nor diagonaliz- 
able by a rotation, nor do its (left or right) eigenvectors define an orthonormal coordinate 
frame (for brevity we omit in what follows the time superscript, M k — > M). As any matrix 
with real elements, M can be expressed in the singular value decomposition form 

M = UDV T , (68) 

where D is diagonal and real, and U, V are orthogonal matrices. The diagonal elements 
<7i, 02, ■ ■ ■ , Od of D are called the singular values of M, namely the square root of the 
eigenvalues of M T M = VD 2 V T (or MM T = UD 2 U T ), which is a symmetric, positive 
semi-definite matrix (and thus admits only real, non-negative eigenvalues). 

Singular values {a,} are not related to the M eigenvalues {Aj} in any simple way. 
From a geometric point of view, when all singular values are non-zero, M maps the 
unit sphere into an ellipsoid, figure 1 1 (b): the singular values are then the lengths of 
the semiaxes of this ellipsoid. Note however that the singular vectors of M T M that 
determine the orientation of the semiaxes are distinct from the M eigenvectors {e^}, 
and that M T M satisfies no semigroup property along the flow. For this reason the M 
eigenvectors {e^} are sometimes called 'covariant' or 'covariant Lyapunov vectors', in 



order to emphasize the distinction between them and the singular value decomposition 
semiaxes directions. 

Eigenvectors / eigenvalues are suited to study of iterated forms of a matrix, such as 
M k or exponentials exp(fA), and are thus a natural tool for study of dynamics. Singular 
vectors are not. They are suited to study of M itself, and the singular value decomposition 
is convenient for numerical work (any matrix, square or rectangular, can be brought to 
this form), as a way of estimating the effective rank of matrix M by neglecting the small 
singular values. 



A.l. Deterministic state space partitions 

We streamline the notation by introducing local coordinate systems z a centered on 
the trajectory points x a , together with a trajectory-centered notation for the map (59), its 
derivative, and, by the chain rule, the derivative (62) of the kth iterate f k evaluated at the 
point x a , 

X = X a -YZa-, fa(Za) = f(x a + Z a ) 
Za+l = M aZa + --- (69) 

M a = f(x a ), M k a =M a+k _ v --M a+l M a , k>2. 
The monodromy (or Floquet) matrix, 

M p , a =M n p(x a ), (70) 

evaluated on the periodic point x a is position dependent, but its eigenvalues, the Floquet 
multipliers {A Pj i, A P; 2> ^p,d} are invariant, intrinsic to the periodic orbit. 

For example, if / is a 1 -dimensional map, its Taylor expansion about x a = x — z a is 

f(x) =X a+l +f a Z a +^faZl+--- » ( 71 ) 

where 

fa = /'(*«), fa=f"W) 

fa = fa+k-Vfa+lfa, k>2. (72) 

A cycle point for which f a = is called a critical point. 

We label trajectory points by either x n , n = 1,2, • • • , in order to emphasize as time 
evolution of xq, or by x a , to emphasize that the trajectory point lies in the state space 
region labeled 'a.' Then the label a+l is a shorthand for the next state space region b 
on the orbit of x a , Xb = x a +i = f(x a ). For example, in figure 11 (a) a periodic point is 
labeled a = Oil by the itinerary with which it visits the regions of the partitioned state 
space ^# = {^#o 5 } i an d as x\ \q = /(xoi l ) > the next point label is b = 1 10. The whole 
periodic orbit Oil = (xon,xno,xioi) ^ s traversed in 3 iterations. 
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FIGURE 11. (a) Periodic points of the map f(x) = 6x(l — x)(l — 0.6x) labeled according to the 
partition j$ = {^#o,^#i}. (b) For a prime cycle p, Floquet matrix M p returns an infinitesimal spherical 
neighborhood of xq £ ^ p stretched into a cigar-shaped ellipsoid, with principal axes given by the eigen- 
directions eW of M p , the monodromy matrix (70). 



A.2. Periodic orbit theory 

Since its initial formulations by Ruelle [67] and Gutzwiller [76], the periodic orbit 
theory has developed into a powerful theoretical and computational tool for prediction 
of quantities measurable in chaotic dynamics. Schematically (a detailed exposition can 
be found in refs. [24, 82]), one of the tasks of a theory of chaotic systems is to predict 
the long-time average of an experimentally measurable quantity a(x) from the spatial 
and time averages 

(a) = lim - (A n ) , A n {xo) = £ a(x k ) . (73) 

What makes evaluation of such averages difficult is chaotic dynamics' sensitivity 
to initial conditions; exponentially unstable trajectories can be tracked accurately only 
for finite times. The densities of trajectories p(x,t), however, can be well behaved for 
t — > °° . Hence the theory is recast in the language of linear evolution operators (Liouville, 
Perron-Frobenius, Ruelle- Araki, • • ■ ) 

p(y,t) = [S% etP ](y)= f dx8b-f(x))p{x,0), (74) 

This evolution operator assembles the density p(y,t) at time t by going back in time to 
the density p (x, 0) at time t = 0. Here we shall refer to the integral operator with singular 
kernel (74) 

& det {*>y) = 8{*-f(y)) 05) 



as the Perron- Frobenius operator, with the subscript det indicating that this is determin- 
istic, to distinguish it from the noisy Fokker-Planck operator (3). 

The Koopman operator action on a state space function a(x) is to replace it by its 
downstream value time t later, a(x) — > a(x(t)) evaluated at the trajectory point x(t): 

[2&a](x) = a(f(x))= f dyJ?tl t (x,y)a(y) 
^dlt^y) = S(y-f(x)). (76) 
Given an initial density of representative points p{x), the average value of a(x) evolves 



as 



= -r^—rf dxdya(y)8(y-f(x)) p(x). 
\PJ(\JJ( 

The 'propagator' 8{y — f (x) ) can equally well be interpreted as belonging to the Perron- 
Frobenius operator (75), so the two operators are adjoint to each other, 

[ dx[^l t a](x)p(x) = I ' dya{y)[^ det p]{y). (77) 

This suggests an alternative point of view, which is to push dynamical effects into the 
density. In contrast to the Koopman operator which advances the trajectory by time 
t, the Perron-Frobenius operator depends on the trajectory point time t in the past. 
Koopman operators are so cool, that it is no wonder that Igor Mezic is so enamored 
with Koopmania [83, 84]. 

The Perron-Frobenius operators are non-normal, not self-adjoint operators, so their 
left and right eigenvectors differ. The right eigenvectors of a Perron-Frobenius oper- 
ator are the left eigenvectors of the Koopman, and vice versa. While one might think 
of a Koopman operator as an 'inverse' of the Perron-Frobenius operator, the notion of 
'adjoint' is the right one, especially in settings where flow is not time-reversible, as is the 
case for infinite dimensional flows contracting forward in time and for stochastic flows. 

If the map is linear, f(x) — Ax, the Perron-Frobenius operator is 



f 1 x 

&det°p{x) = J dyd(x-Ay)p(y) = ^Tp(r^r) 



(78) 



In the |A| > 1 expanding case the right, expanding deterministic eigenfunctions [12, 85] 
are monomials 

p k {x)^x k /k\, * = 0,1,2,-.., (79) 
with eigenvalues 1/|A|A", while the left, contracting eigenfunctions are distributions 



Pi W^(-lf5«(i). 



(80) 



In discretizations J£\ et (y,x) is represented by a matrix with y, x replaced by dis- 
crete indices, and integrals over x replaced by index summation in matrix multiplica- 
tion. Indeed, for piece-wise linear mappings Perron-Frobenius operator can be a finite- 
dimensional matrix. For example, consider the expanding 1 -dimensional 2-branch map 
f(x) with slopes Ao > 1 and Ai = — Ao/(Ao — 1) < — 1: 

f ( r \- f M x ) = A ox> xeJZ = [0,1/Ao) S1 

/i(jc) = Ai(1-x), xe^! = (l/Ao,l]. 

As in figure 1 1 (a), the state space (i.e., the unit interval) is partitioned into two regions 
= {^#o,^i }• If density p(x) is a piecewise constant on each partition 



the Perron-Frobenius operator acts as a [2 x 2] Markov matrix L with matrix elements 

*)->Lp=( ^ (83) 

stretching both po and pi over the whole unit interval A. Ulam [47, 2] had conjectured 
that successive refinements of such piece-wise linear coarse-grainings would provide 
a convergent sequence of finite-state Markov approximations to the Perron-Frobenius 
operator. 

The key idea of the periodic orbit theory is to abandon explicit construction of natural 
measure (the density functions typically observed in chaotic systems are highly singular) 
and instead compute chaotic spatial and time average (73) from the leading eigenvalue 
zo = z(j3) of an evolution operator by means of the classical trace formula [86, 24], 
which for map /, takes form 



{a) =dj5 



1 - z ™ Pe rP-A p 

I^I iwi-^r (84) 



/5=o ' ^o z ~ z <x p r=i|det(l-M; 
or, even better, by deploying the associated spectral determinant 

(°° 1 7 rn P e r P- A P 

These formulas replace the chaotic, long-time uncontrollable flow by its periodic orbit 
skeleton, decomposing the dynamical state space into regions, with each region M p 
centered on an unstable periodic orbit p of period n p , and the size of the p neighborhood 
determined by the linearization of the flow around the periodic orbit. Here M p is the 
monodromy matrix (62), evaluated in the periodic orbit p, the deterministic exponential 
contraction/expansion is characterized by its Floquet multipliers {A^i,--- ,A p j}, and 
p contribution to (84) is inversely proportional to its exponentiated return time (cycle 
period n p ), and to the product of expanding eigenvalues of M p . With emphasis on 



expanding: in applications to dissipative systems such as fluid flows there will be only 
several of these, with the contracting directions - even when their number is large or 
infinite - playing only a secondary role. 

Periodic solutions (or 'cycles') are important because they form the skeleton of the 
invariant set of the long time dynamics [9, 87], with cycles ordered hierarchically; short 
cycles give dominant contributions to (84), longer cycles corrections. Errors due to ne- 
glecting long cycles can be bounded, and for hyperbolic systems they fall off expo- 
nentially or even super-exponentially with the cutoff cycle length [31]. Short cycles 
can be accurately determined and global averages (such as transport coefficients and 
Lyapunov exponents) can be computed from short cycles by means of cycle expan- 
sions [9, 87, 82,24]. 

A handful of very special, completely hyperbolic flows are today mathematically fully 
and rigorously under control. Unfortunately, very few physically interesting systems are 
of that type, and the full picture is more sophisticated than the cartoon (84). 



B. FOKKER-PLANCK OPERATOR, CONTINUOUS TIME 

FORMULATION 

The material reviewed in this appendix is standard [1, 3, 14], but needed in order to set 
the notation for what is new here, the role that Fokker-Planck operators play in defining 
stochastic neighborhoods of periodic orbits. 
Consider a J-dimensional stochastic flow 

f = v(*) + |(f), (86) 

where the deterministic velocity field v(x) is called 'drift' in the stochastic literature, 
and q(t) is additive noise, uncorrelated in time. A way to make sense of q(t) is to 
first construct the corresponding probability distribution for additive noise t, at a short 
but finite time 8t. In time 8t the deterministic trajectory advances by v(x n ) 8t. As 8t 
is arbitrary, it is desirable that the diffusing cloud of noisy trajectories is given by a 
distribution that keeps its form as 8t — > 0. This holds if the noise is Brownian, i.e., the 
probability that the trajectory reaches x n+ \ is given by a normalized Gaussian 



28t^"A^ 



(87) 



Here % n — 8x n — v(x n ) 5?, the deviation of the noisy trajectory from the deterministic 
one, can be viewed either in terms of velocities {x,v(x)} (continuous time formulation), 
or finite time maps {x n — > x n+ i,x n — > f 5t (x n )} (discrete time formulation), 



8x n = x n+i -x n ~x n 8t, f 5t (x n )-x n ~v(x n )8t , (88) 



where 

{x ,xi,--- ,*„,-•• i x k} = {x(0),x(5t),--- ,x(n8t),--- ,x(t)} (89) 



is a sequence of k+ 1 points x n = x(t n ) along the noisy trajectory, separated by time in- 
crements 8t = t/k, and the superfix T indicates a transpose. The probability distribution 
t, (t n ) is characterized by zero mean and covariance matrix (diffusion tensor) 

(Zj(t n )) = 0, (S i (t m )Sj(t n )) = A iJ 8 nm , (90) 

where (■ ■ ■) stands for ensemble average over many realizations of the noise. For exam- 
ple, in one dimension the white noise t, n = x n+ i —x n for a pure diffusion process (no 
advection, v(x n ) = 0) is a normally distributed random variable, with standard normal 
(Gaussian) probability distribution function, 



^ (j: ' xo)= vib exp 



(x-XqY 



(91) 



2At 

of mean 0, variance At, and standard deviation \/At, uncorrected in time: 

(x n+ i -x n ) =0, (Cwi -x m )(x n+ i -x n )) =A8 mn . (92) 

Jz?'(x,xo) describes the diffusion at any time, including the integer time increments 
{t„} = {St, 2St, ■ ■ ■ ,nSt, ■■■}, and thus provides a bridge between the continuous and 
discrete time formulations of noisy evolution. We have set 8t = 1 in (92) anticipating 
the discrete time formulation of sect. 2. 

In physical problems the diffusion tensor A is almost always anisotropic: for example, 
the original Langevin flow [88] is a continuous time flow in {configuration, velocity} 
phase space, with white noise probability distribution exp(— x 2 /2kBT) modeling random 
Brownian force kicks applied only to the velocity variables x. In this case one thinks 
of diffusion coefficient D = kgT/2 as temperature. For sake of simplicity we shall 
sometimes assume that diffusion in d dimensions is uniform and isotropic, A(x) = 2D1. 
The more general case of a tensor A which is a state space position dependent but time 
independent can be treated along the same lines, as we do in (14). In this case the 
stochastic flow (86) is written as [89] dx = v(x) dt + o(x) dt, (t) , o{x) is called 'diffusion 
matrix,' and the noise is referred to as 'multiplicative.' 

The distribution (87) describes how an initial density of particles concentrated in a 
Dirac delta function dXx n spreads in time St. In the Fokker- Planck description individual 
noisy trajectories are replaced by the evolution of the density of noisy trajectories. The 
finite time Fokker-Planck evolution p(x,t) = J??' op(j,0) of an initial density p(^o,0) 
is obtained by a sequence of consecutive short-time steps (87) 



^(x^xq) = J [dx] exp | - £ [x n+l -f 5t (x n )} 2 I , (93) 



where t = k8t, and the Gaussian normalization factor in (87) is absorbed into interme- 
diate integrations by defining as 

k-l 

[dx] = N- l Y[dxi 

n=l 

N = (2^5?)^/ 2 (detA) 1//2 anisotropic diffusion tensor A 

= {27tA8t) d/2 isotropic diffusion, (94) 

The stochastic flow (86) can now be understood as the continuous time, 8t — > limit, 
with the velocity noise £, (t) a Gaussian random variable of zero mean and covariance 
matrix 

%(t)) =0, (&t)Zj(t')) =Aij8(t-t'). (95) 



It is worth noting that the continuous time flow noise f (t) in (86) and (95) is dimension- 
ally a velocity [*]/[*], while the discrete time noise E, n in (87), (90) is dimensionally a 
length [x]. The continuous time limit of (93), 8t = t/k — > 0, defines formally the Fokker- 
Planck operator 

J? t (x,x ) = J[dx]exp^-^J\x(x)-v(x(t))] 2 dx^ (96) 

as a stochastic path (or Wiener) integral [73, 90, 3] for a noisy flow, and the associated 
continuous time Fokker- Planck (or forward Kolmogorov) equation [1, 3, 91] describes 
the time evolution of a density of noisy trajectories (86), 

dtp(x,t) + V-(y(x)p(x,t)) = DV 2 p{x,t). (97) 

The 8t — > limit and the proper definition of are delicate issues [92, 93, 14, 94] of 
no import for the applications of stochasticity studied here. 

In probabilist literature [95] the differential operator —V • (v(x)p(x,t)) + DV 2 p(x,t) 
is called 'Fokker- Planck operator;' here we reserve the term exclusively for the finite 
time, 'Green function' integral operator (96). The exponent 



1 



2A8t 



x n+1 -f St (x n )] 2 ~ --L[i(T)-v(x(T))] 2 8t (98) 



can be interpreted as a cost function which penalizes deviation of the noisy trajectory 
8x from its deterministic prediction v8t, or, in the continuous time limit, the deviation 
of the noisy trajectory tangent x from the deterministic velocity v. Its minimization is 
one of the most important tools of the optimal control theory [96, 97], with velocity x(x) 
along a trial path varied with aim of minimizing its distance to the target v(x(x)). 

The finite time step formulation (93) of the Fokker-Planck operator motivates the 
exposition of sect. 2, which starts by setting 8t = 1. In the linearized setting, the two 
formulations are fully equivalent. 



Remark B.l A brief history Of noise. The cost function (98) appears to have been first 
introduced by Wiener as the exact solution for a purely diffusive Wiener-Levy process in one 



dimension, see (91). Onsager and Machlup [98, 99] use it in their variational principle to 
study thermodynamic fluctuations in a neighborhood of single, linearly attractive equilibrium 
point (i.e., without any dynamics). The dynamical 'action' Lagrangian in the exponent of (96), 
and the associated symplectic Hamiltonian were first written down in 1970's by Freidlin and 
Wentzell [99], whose formulation of the 'large deviation principle' was inspired by the Feynman 
quantum path integral [100]. Feynman, in turn, followed Dirac [101] who was the first to 
discover that in the short-time limit the quantum propagator (imaginary time, quantum sibling of 
the Wiener stochastic distribution (91)) is exact. Gaspard [4] thus refers to the 'pseudo-energy of 
the Onsager-Machlup-Freidlin-Wentzell scheme.' M. Roncadelli [22, 102] refers to the Fokker- 
Planck exponent in (96) as the 'Wiener-Onsager-Machlup Lagrangian,' constructs weak noise 
saddle-point expansion and writes transport equations for the higher order coefficients. 



B.l. Ornstein-Uhlenbeck process 

The variance (17) is stationary under the action of Jzf, and the corresponding Gaussian 
is thus an eigenfunction. Indeed, for the linearized flow the entire eigenspectrum is 
available analytically, and as Q a can always be brought to a diagonal, factorized form in 
its orthogonal frame, it suffices to understand the simplest case, the Ornstein-Uhlenbeck 
process in one dimension. This simple example will enable us to show that the noisy 
measure along unstable directions is described by the eigenfunctions of the adjoint 
Fokker-Planck operator. 

The simplest example of a stochastic flow (86) is the Langevin flow in one dimension, 

^=Ax + |(0, (99) 

with 'drift' v(x) linear in x, and the single deterministic equilibrium solution x = 0. 
The associated Fokker-Planck equation (97) is known as the Ornstein-Uhlenbeck pro- 
cess [68, 11, 12, 103, 3]: 

d t p(x,t) + d x (Xxp(x,t)) =Dd?p(x,t). (100) 

(Here A = 2D, D = Einstein diffusion constant.) One can think of this equation as 
the linearization of the Fokker-Planck equation (97) around an equilibrium point. For 
negative constant A the spreading of noisy trajectories by random kicks is balanced by 
the linear damping term (linear drift) v(x) = Xx which contracts them toward zero. For 
this choice of v(x), and this choice only, the Fokker-Planck equation can be rewritten 
as the Schrodinger equation for the quantum harmonic oscillator, with its well-known 
Hermite polynomial eigenfunctions [104, 105] discussed here in sect. 4.1. 

The key ideas are easier to illustrate by the noisy, strictly equivalent discrete-time 
dynamics of sect. 4.1, rather than by pondering the meaning of the stochastic differential 
equation (100). 

Remark B.2 Quantum mechanical analogue. The relation between Ornstein-Uhlenbeck 
process and the Schrodinger equation for the quantum harmonic oscillator is much older than 



quantum mechanics: Laplace [13] wrote down in 1810 what is now known as the Fokker- 
Planck equation and computed the Ornstein-Uhlenbeck process eigenfunctions [106] in terms 
of Hermite polynomials (33). According to L Arnold [14] review of the original literature, the 
derivations are much more delicate: the noise is colored rather than Dirac delta function in time. 
He refers only to the linear case (99) as the 'Langevin equation'. 



C. LYAPUNOV EQUATION 

In his 1892 doctoral dissertation [26] A. M. Lyapunov defined a dynamical flow to be 
"stable in the sense of Lyapunov" if the scalar Lyapunov function 

V(x)>0, V(0) = 

computed on the state space of dynamical flow x = v(x) satisfies the 'inflow' condition 

dV 



V 



dx~j Vj -°' (101) 



While there is no general method for constructing Lyapunov functions, for a linear d- 
dimensional autonomous flow, x = Ax, the Lyapunov function can be taken quadratic, 



and (101) takes form 



V(x)=x T ^x, Q = Q T >0, (102) 



V x' Ia' ' • 1 ,\ I.Y. 



Q Q 



Here A T is the transpose of the stability matrix A, and Q > is the shorthand for matrix 
being positive definite (or Hurwitz), i.e., having the entire eigenvalue spectrum, o(Q) G 
C+, in the right-hand half of the complex plane. Strict positivity Q > guarantees that 
Q is invertible. 

The Lyapunov differential equation for a time varying system is 

Q = AQ + QA T +A, Q(t ) = Q . (103) 

For steady state, time invariant solutions Q = 0, (103) becomes the Lyapunov matrix 
equation 

AQ + QA T +A = 0. (104) 

The Lyapunov theorem [26, 108] states that a [dxd] matrix Q has all its characteristic 
roots with real parts positive if, and only if, for any positive definite symmetric matrix 
A = A T > 0, there exists a unique positive definite symmetric matrix Q that satisfies the 
continuous Lyapunov equation (104). The flow is then said to be asymptotically stable. 
In our application we are given a noise correlation matrix A, and the theorem states 
the obvious: the stationary state with a covariance matrix Q > exists provided that the 



stability matrix A is stable, A < 0. We have to require strict stability, as A < would allow 
for a noncompact, Brownian diffusion along the marginal stability eigen-directions. 

For a linear discrete-time system x n+ \ = Mx n the quadratic Lyapunov function (102) 
must satisfy 

V(x n+1 ) -V(x n ) =x T n (m t ^M- ^ x n < 0, 

leading to the discrete Lyapunov equation 

Q = MQM T + A , for any A = A T > 0. (105) 

Lyapunov theorem now states that there is a unique solution Q, provided the Jacobian 
matrix M is a convergent matrix, i.e., a matrix whose eigenvalues (Floquet multipliers) 
are all less than unity in magnitude, |A,-| < 1 . 

We note in passing that the effective diffusive width is easily recast from the discrete 
map formulation back into the infinitesimal time step form of appendix B: 

M = e A8 \ M T = e AT§ \ A — >■ <5/A, (106) 

where A = dv/dx is the stability matrix. Expanding to linear order yields the differential 
version of the equilibrium condition (17) 

0=AQ + QA T +A (107) 

(with the proviso that now A is covariance matrix (95) for the velocity fluctuations). 
The condition (107) is well known [14, 109, 3] and widely used, for example, in the 
molecular and gene networks literature [110, 111, 112, 113]. In one dimension, A — > A 
(see (63)), A diffusion tensor — > A, and the diffusive width is given by [14] 

Q = -A/2X, (108) 

a balance between the diffusive spreading A and the deterministic contraction rate X < 0. 
If A — > 0, the measure spreads out diffusively. 

If A has eigenvalues of both signs, the necessary and sufficient condition for the 
existence of a unique solution to time invariant case (104) is that no two eigenvalues 
add up to zero [114], 

A^+A^VO, i,j= 1,2,- •• ,d. (109) 

For the discrete Lyapunov equation (105) the condition for the existence of a unique 
solution is that no two eigenvalues have product equal to one, 

AfAy^l, i,j=l,2,---,d. (110) 

This condition is obviously satisfied in the Lyapunov case, with Jacobian matrix M 
asymptotically stable in the discrete-time domain (all eigenvalues of M are strictly inside 
of a unit circle). 



If A is stable, the continuous Lyapunov equation (104) has a unique solution 



poo 

Q= / dte At Ae A \ (111) 
Jo 

and 

oo 

Q=Y,M k A(M T ) k 

k=0 

is the unique solution of the discrete Lyapunov equation (105). 

In chaotic dynamics we are interested in saddles, i.e., hyperbolic points with both 
expanding and contracting eigen-directions. Given a [dxd] square matrix A with real 
elements, let the numbers of eigenvalues in the left half of the complex plane, the 
imaginary axis and the right half of the complex plane be denoted by the integer triple 

In A = (InA_,InAo,InA + ) 

which counts (stable, marginal, unstable) eigen-directions. Following J. J. Sylvester 
(1852), this triple is called the inertia of matrix A. In the case of a nondegenerate 
symmetric bilinear form (such as the symmetric matrix Q) the numbers of positive 
and negative eigenvalues are also known as the signature of the matrix. A is said to 
be (negative) stable if A < 0, i.e., InA = (J, 0,0), and positive stable if A > 0, i.e., 
InA = (0,0,J). 

The Lyapunov theorem generalized to hyperbolic fixed points is known as the Main 
Inertia Theorem [115, 116, 117]: For a given stability matrix A and a noise correlation 
matrix A = A T > 0, there exists a symmetric Q such that AQ + QA T + A = and 

(In(2_,0,In(2 + ) = (InA + ,0,InA_), 

if and only if A has no marginal eigenvalues [118], InAo = 0. The Lyapunov theorem is 
a special case: it states that Q > if and only if InA = (d, 0, 0) . 

We shall solve (104) and diagonalize Q numerically. The Main Inertia Theorem 
guarantees that the covariance matrix Q will have the same number of expanding / 
contracting semi-axes as the number of the expanding / contracting eigenvalues of the 
stability matrix A. The contracting ones define the semi-axes for covariance evolution 
forward in time, and the expanding ones the semi-axes for the adjoint evolution. 

Remark C.l Lyapunov equation. The continuous Lyapunov equation (104) is a special 
case of the Sylvester equation, 

AQ + QB = C (112) 

where A, B, Q, C are [dxd] matrices, and the discrete Lyapunov equation (105) is a special case 
of Stein's equation. The Sylvester equation can be solved numerically with the Bartels-Stewart 
algorithm [119], in full generality, with no assumptions on degenerate eigenvalues or defective 
matrices (matrices for which there are fewer eigenvectors than dimensions). It is implemented 
in LAPACK, Matlab and GNU Octave. Kuehn [18] reviews the available numerical meth- 
ods for solving the Lyapunov equation. Discrete Lyapunov equation is solved numerically with 
the Kitagawa algorithm [120], using, for example, Mathematica DiscreteLyapunovSolve. 



The Mathematica solvers for the continuous- and discrete-time Lyapunov equation implement 
the Schur method by Bartels and Stewart [119], based on the decomposition of to the real 
Schur form; and the Hessenberg-Schur method by Golub, Nash, and Van Loan [121] to solve 
the Sylvester equation, based on the decomposition of the smaller of two matrices and to the 
real Schur form and the other matrix to the Hessenberg form [122]. For continuous Lyapunov 
equation see Matlab function lyap, which cites ref. [123, 124, 125, 126, 127, 121] as sources. 
Related is Matlab function covar for output and state covariance of a system driven by white 
noise. For a more practical criterion that matrix has roots with negative real parts, see Routh- 
Hurwitz stability criterion [128] on the coefficients of the characteristic polynomial. The crite- 
rion provides a necessary condition that a fixed point is stable, and determines the numbers of 
stable/unstable eigenvalues of a fixed point. 



D. CONFLUENT HYPERGEOMETRIC FUNCTIONS AND 
LAGUERRE POLYNOMIALS 

Let now transform this new density, around the next pre-image x a -i = f~ l {x a -\), 



JS?V°^-i = e a y\dy] (113) 

where a 2 = f" l _ 1 2 /S(Q a + A). Now change the variable £ = y^faJlA, and write the 
density p a -\ (y) as a power series, so that the previous integral reads 



= M / W1 r(*-^)' fur [(2A)yr 
5— <^*-*> 1(4^3^ U/^- 2 r-J <114) 

We then group all the terms up to order 0(A) and neglect 0(A 2 ) and higher, and call 

V = Vaf' a _ 2 z a -2 

h " ! t Q 4[n\(4n-2)\] 1 

e^ 4 -2A[3ar ? 2 ]4»Q,^_T 7 ^ (115) 

where 4> (also sometimes called M or 1F1 in the literature) is a confluent hypergeometric 
function of the first kind [107], which can be expressed as 



We now want to evaluate the variance of the density (115) 

J dZa-2Z 2 a _ 2 pa-2(Za- 2 ) 

Qa-2 = Ya ~ 1 \~ ' " ll > 

J dZa-2Pa-2{Za-2) 

It is useful to know, when computing the denominator of (117), that 

J7 ] 2 ^{l 3 -^d7]=0 (118) 

so that 

(« 2 ^ 4 _ 2 )- 1 / 2 jT > 2 e -^T > -4D[3/ fl l|]jT > 4 4>(| > | -t; 4 )^ 

- 1 r (3/4) 1 I A)- Qa ~ l+A (119) 

in the last identity we used the definition of a and (50). 

Next we derive (116), which expresses the a confluent hypergeometric function in 
terms of an exponential and a Laguerre polynomial. Start with the identity [107]: 

<f>(a,b,z) =e z <Z>(b-a,b,-z) (120) 
in particular, the hypergeometric function in (1 15) becomes 

*g.i.-n 4 )=.-'**(-i,^ 4 ) (12D 

A confluent hypergeometric function can be written in terms of a Laguerre polynomial 
[107]: 

L«(x)=f" + a )^(-n,a + l,x),, L«(x) = £ (-irf" + a )^ (122) 
n / ^ \n-mj ml 



Thus, in our case 



and 

* (-I, |,i| 4 ) *-" 4 = 1 C| - 1|- ) . (124) 
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